Self- similar Approximants of the Permeability 
in Heterogeneous Porous Media from Moment Equation Expansions 



S. Gluzman 1 and D. Sornette 1 ' 2 ' 3 
1 Institute of Geophysics and Planetary Physics 
University of California Los Angeles, Los Angeles, CA 90095-1567 

2 Department of Earth and Space Sciences, UCLA 

3 Laboratoire de Physique de la Matiere Condensee 

CNRS UMR 6622 and Universite de Nice-Sophia Antipolis, 06108 Nice Cedex 2, France 

February 1, 2008 



Abstract 

We use a mathematical technique, the self-similar functional renormalization, to construct formulas for the 
average conductivity that apply for large heterogeneity, based on perturbative expansions in powers of a small pa- 
rameter, usually the log-variance a\- of the local conductivity. Using perturbation expansions up to third order and 
fourth order in ay obtained from the moment equation approach, we construct the general functional dependence of 
the transport variables in the regime where o\ is of order 1 and larger than 1. Comparison with available numerical 
simulations give encouraging results and show that the proposed method provides significant improvements over 
available expansions. 
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1 Introduction 



Subsurface hydraulic parameters such as medium permeability and porosity, saturation curves and relative permeabil- 
ity have been traditionally viewed as well-defined local quantities that can be assigned unique values at each point in 
space. Yet, subsurface flow takes place in a complex environment whose makeup varies in a manner that cannot be 
predicted deterministicaily in all of its relevant details. This makeup tends to exhibit discrete and continuous varia- 
tions on a multiplicity of scales, causing hydraulic parameters to do likewise. In practice, such parameters can at best 
be measured at selected locations and depth intervals where their values depend on the scale (support volume) and 
mode (instrumentation and procedure) of measurement. Estimating the parameters at points where measurements are 
not available entails a random error. Quite often, the support of measurement is uncertain and the data are corrupted 
by experimental and interpretive errors. 

Uncertainty is usually dealt with either deterministicaily through upscaling or stochastically through the evaluating 
statistical moments. Statistical moments can be obtained through Monte Carlo simulations or development of moment 
differential equations, which is the method from which we start our analysis. 

In the stochastic approach, parameter values determined at various points within a more-or-less distinct soil unit 
can be viewed as a sample from a random field defined over a continuum. This random field is characterized by a 
joint (multivariate) probability density function or, equivalently, its joint ensemble moments. Thus, a parameter such 
as (saturated, natural) log hydraulic conductivity Y(x) = hiK s (x) varies not only across the real space coordinates 
x within the unit, but also in probability space (this variation may be represented by another "coordinate" £, the 
configuration coordinate, which, for simplicity, we suppress). Whereas spatial moments are obtained by sampling 
y(x) in real space (across x), ensemble moments are defined in terms of samples collected in probability space 
(across £). 

In the moment equation approach, which we propose to exploit here, the stochastic differential equations are 
averaged first to obtain moments differential equations (MDEs) governing the statistical moments of the dependent 
variables. The MDEs are themselves deterministic and can be solved numerically or sometimes, analytically. The 
MDE approach has important advantages. First, only a small number of equations must be solved: one for the mean 
and one each for a small number of variances and covariances. Second, the coefficients of the MDEs are relatively 
smooth because they are averaged quantities. Thus the MDEs can be solved on comparatively smooth grids. Third, the 
MDEs are available in analytical form, even though they are usually solved numerically in applications. This holds the 
potential for increased physical understanding of the mechanisms of uncertainty through qualitative analysis. Finally, 
in many applications MDE approaches provide a good estimate of the behavior of large variance systems despite 
being based on small perturbation theory. 

Since the moment equation approach derives the equations of evolutions of the moments of the distribution of 
the transport variables by averaging the stochastic differential equations of transport in heterogeneous porous me- 
dia, its fundamental limitation is the assumption that the variance a\ of the log hydraulic conductivity is small, 
in contradiction with real geological settings of interest to a large variety of geophysical applications. In order to 
obtain reliable descriptions of the transport coefficients for large <7 y , our goal in the present paper is to adapt and ex- 
tend the self-similar functional renormalization method to the moment equation approach. In essence fundamentally 
non-perturbative, this recently developed technique provides us with a stable and robust estimation of the transport 
variables at large values of the perturbation parameter a\. The functional renormalization method associates ideas 
from the renormalization group theory of multiscale and critical phenomena [21] with methods from the theory of 
dynamical systems and of control theory. Using perturbation expansions up to third order and fourth order in <7y ob- 
tained from the moment equation approach, we construct the general functional dependence of the transport variables 
in the regime where a\ is of order 1 and larger than 1. 

The next section 2 recalls briefly how perturbation expansions are obtained from the moment equation approach. 
Section 3 summarizes the general formulation of the self-similar approximation theory. Section 4 gives the results of 
the application of the functional renormalization method to the moment equation expansions at increasing orders in 
Oy. Section 5 formulates the expansion in powers of 1/d, where d is the space dimension. This procedure well-known 
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in statistical physics is resummed by the fonctional renormalization to provide accurate formulas. Section 6 briefly 
outlines future directions of investigations. 



2 Problem Formulation and perturbation expansions 
2.1 Basic equations 

It has become common to quantify uncertainty in ground water flow models by treating hydraulic conductivity, K, 
and derived quantities like hydraulic head, h, as random fields. For steady-state flows in the absence of sources and 
sinks, the statistics of h can be obtained from the stochastic flow equation 

V • [K (x) V/i(x)] = (1) 

when the statistics of K are known. We further assume that the site of interest is sufficiently characterized so that 
available experimental data are sufficient to obtain the statistics of K, such as its (ensemble) mean, K, variance, a\, 
and (two-point) correlation structure, px(x, y). Then one can solve directly for the moments of h by developing 
deterministic equations for the moments from ( [[]). In general this involves taking the expected value of ([[]) and 
similar equations for higher-order moments, closing the system of moment equations (usually through perturbation 
approximations). Numerical solutions for moment equations are typically computationally more efficient than Monte 
Carlo simulations. In the first place, taking expected values smoothes parameters in the moment equations which 
in turn allows low-resolution grids for numerical solutions. Furthermore, the number of moment equations is much 
smaller than the number of realizations required by Monte Carlo simulations. Additionally, the moment equations 
lend themselves to qualitative analysis. 

Here, we concentrate on flow through highly heterogeneous porous media with the variance a\ of log hydraulic 
conductivity Y = \riK as the expansion parameter of the theory. We estimate the mean hydraulic head, /i(x), and 
assess the errors associated with such an estimation. We represent K(x) = K(x) + K'(x.) as the sum of a mean, 
K(x), and a zero-mean random deviation, K'(x), with variance a\ (x). Similarly, /i(x) = /i(x) + h'(x) with 
h'(x) = and variance cr?(x). 

The average steady-state flow equation becomes 

V • [f(x)V%)] + V • r(x) = (2) 

which consists of a deterministic mean part, KVh, and a deterministic residual flux, r = —K'Vh'. Solutions of (Q) 
require the mean conductivity, K (x), and in most cases, a method for closing an expansion of r(x). Usually r(x) is 
approximated through perturbation expansions based on ay, the variance of Y = In K, the logarithm of conductivity. 
This approach works well as long as a\ is small. This restriction is a stumbling block on the road to applicability 
of numerous theoretical analyses to real-world problems. Our goal here is to provide a general theoretical method to 
extend the domain of application of moment equations to the large heterogeneity limit. 



2.2 Perturbation Expansions 

Consider asymptotic expansions of the parameters and functions, K = K g (1 + a\ /2 + ...); q = + q^ + • . . 

; h = + h + . . . ; and r = + . . ., where K g = exp(Y), Y being the ensemble mean of Y. The superscript 
(i) denotes terms that are of ith-order, i.e. contain only the ith power of a\. The first-order (in Oy) approximation of 
the residual flux is given by [22, and references therein] 

•W(x) = K g o\\ /3 y(y,x)V x V^G(y,x)V^ ( ° ) (y)(iy, (3) 



where py(y, x) is the spatial two-point autocorrelation function of Y, and G(y, x) is the deterministic Green's func- 
tion for Laplace equation in Q, subject to the corresponding homogeneous boundary conditions. It is a standard 
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practice in stochastic hydrogeology to rely on the first-order approximation of r [0J, but higher-order approximations 
are also available [13]. 

Collecting the terms of the same powers of a\ yields the zeroth-order approximation of the mean head in |2|, 

^ g V 2 ^ (0) (x) =0, (4) 

and its first-order approximation, 

m Tit?, —ir\\ 

= 0. (5) 



KgV 2 h {l) (x) + V 



K 5 V/i w (x) - r W(x) 



Solving a system of these sequential approximations leads to h = h + h . Strictly speaking, for such expan- 
sions to be asymptotic it is necessary that a\ <C 1, i.e. that porous media be mildly heterogeneous. However, various 
numerical simulations (e.g., [12]) have demonstrated that these first-order approximations remain remarkably robust 
even for strongly heterogeneous media with a\ as large as 4. 



2.3 Effective Conductivity of Porous Media 



For an effective conductivity to exist in the strict sense, it is necessary that V/i be constant. A somewhat a less 
restrictive assumption requires V/i to vary slowly in space, i.e. to have negligibly small derivatives [Qj. Then one can 
localize expression (0) as 



r«(x) « *r fl 4A«(x)V/i (0) (x), A«(x) = f p y (y,x) V X V^ G(y, x) dy 

Jn 



(6) 



Under these conditions, retaining the two leading terms in the asymptotic expansion of the mean Darcy flux, 



q 



[i] 



7(0) 



+ qW, yields 

_ q [1] (x) 
K„ 



V/i (1) (x) + 



I + af (-1 



A«(> 



v7^ (0) (x) 



(7) 



For flow through infinite, statistically homogeneous porous media under mean uniform flow conditions, or at points 

"""" - const and V/T = 



away from boundaries and singularities, the mean hydraulic head gradient J = Vh 
(i > 1) [||> ^3[|- This gives rise to the effective conductivity given approximately by 



K 



ef 



K n 



(8) 



where d is the space dimension. 

Various attempts to generalize this asymptotic expansion to highly heterogeneous formations were attempted by 
conjecturing that expression ([|) represents the two leading terms in the expansion of an exponent [|l5|, [l9| ], 



K, 



ef 



K„ exp 



1 1 

2 ~ d 



(9) 



In recent years the question of validity of expression (g) was the focus of a thorough investigation. It was proven that 
expression @ is rigorously valid under one-dimensional flow in log-normal fields where it yields the harmonic mean 
K h = Kg exp(— o"y/2) [^ [l8[ ]. It is also rigorously valid under two-dimensional flow in log-normal, statistically 
isotropic conductivity fields where it yields the geometric mean K g [15]. For three-dimensional flow in log-normal, 
statistically isotropic fields, the second-order (in o\) term in <js|> was found to be in agreement with the Taylor series 
expansion of @ [|3]] . While unsuccessful! attempts to prove @ for three-dimensional flows in such fields have been 
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reported [14, 17], De Wit [0] demonstrated that the third-order correction in (||) is not equal to the third-order term in 
the Taylor expansion of ( ||), thereby proving this conjecture to be not strictly valid for three-dimensional Gaussian 
isotropic media. Instead, it was demonstrated that this and higher-order terms depend on the shape of the correlation 
function py. 

These results suggest that it would be beneficial to view the equation <JSj> and its higher order terms as a perturbation 
expansion of the true transport variable in powers of the variance a Y . In this sense, the passage from ([|) to @ is a 
resummation procedure. It thus makes full sense to ask what could be the most general and robust resummation that 
can generalize (|8]) in order to extend its domain of validity in the regime of large a Y where the initial perturbation 
expansion breaks down. 

It is often the case that perturbation expansions are not converging but are instead diverging series. Even if the 
series is convergent for small perturbation parameters a Y , one is in general interested in the regime where a Y is 
of order 1 and larger. In this case, the perturbation series is divergent and is of no direct use. The study of such 
summation of divergent series is the problem of great importance in theoretical physics, applied mathematics and 
engineering. This is because realistic problems are usually solved by means of some calculational algorithm often 
resulting in divergent sequence of approximations. Assigning a finite value to the limit of a divergent sequence 
is called renormalization or summation technique. The most widely used such technique is Pade summation [j|]. 
However, the Pade summation method has several shortcomings. First of all, to reach a reasonable accuracy of Pade 
approximants, one needs to possess tens of terms of a perturbation series. In contrast, only a few terms are often 
available because of the complexity of the problem. Second, Pade approximants are defined for the series of integer 
powers. But in many cases asymptotic series arise having noninteger powers. Third, there are quite simple examples 
that are not Pade summable even for a sufficiently small variable. Last but not least, Pade summation is more of a 
numerical technique providing answer in the form of numbers. Therefore, it is difficult, if possible, to analyze the 
results when the considered problem contains several parameters to be varied, since for each given set of parameters 
one has to repeat the whole procedure of constructing a table of Pade approximants and of selecting from them one 
corresponding to a visible saturation of numerical values. 

We thus turn to the method of so-called self-similar approximation or functional renormalization that provide a 
very interesting alternative. We first summarize the idea of the technique and then apply it to calculate properties of 
transport in porous media in the limit of large heterogeneity. 



3 General formulation of the self-similar approximation theory 

General ideas and the mathematical foundation of the self-similar approximation theory have been described in detail 
in [Q, Q 25, 26, [k], [H], ||[ ^g, 29]. The approach is applicable in all cases, when either just a few terms of a 

series are known or when a number of such terms are available. We are always able to obtain analytical formulas that 
are easy to consider with respect to varying characteristic parameters. We now expose the general idea of the method 
of self-similar approximation. 

Consider the case, when for a sought function f(x), one derives an approximate perturbative expansion 

k 

Pk{x) = J^Qn X° n , (10) 
n=0 

in which a n is an arbitrary real number, integer or noninteger, positive or negative. Following the method of the 
algebraic self-similar renormalization ]^J, we define the algebraic transform 

k 

P k (x,s) = x s Pk (x) = ^2a n x s+an , (11) 

n=0 

where s is real. Rather than constructing a trajectory in the functional space of the initial approximations, the idea 
behind the introduction of the transform (O) is to deform smoothly the initial functional space of the approximations 
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Pk{x) in order to obtain a faster and better controlled convergence in the space of the modified functions P k (x, s). 
This convergence can then be mapped back to get the relevant estimations and predictions. The exponent s depends 
on x in general and will be acted upon as a control function in order to accelerate convergence. 

Then, by means of the equation Pq(x,s) = ao x s+a ° = <p, we obtain the expansion function x(ip, s) = 

l/(s+Q!o) _ 

. Substituting the latter into (|10|), we have 



Vk(<P, s) = Pk{x((f, s), s) = S~] a n ( — ) 



(s+a n )/(s+ao) 



(12) 



The family {y k } of transforms (|10|) is called the approximation cascade, since its trajectory {y k ((p, s) \ k = 0, 1, 2...} 
is bijective to the sequence {P k (x,s) [ k = 0, 1,2...} of approximations ( A cascade is a dynamical system 
in discrete time k = 0, 1, 2..., whose trajectory points satisfy the semigroup property y k + p (ip, s) = ykiVpip, s), s). 
The physical meaning of the above semigroup relation can be understood as the property of functional self-similarity 
with respect to the varying approximation number. The self-similarity relation is a necessary condition for the fastest 
convergence criterion. 



For the approximation cascade {yk}, defined by transform (|12|), the cascade velocity is 

(s+a k )/(s+ao) 



v k (ip, s) = y k (ip, s) - Vk-ii'P, s) = a k 
This is to be substituted into the evolution integral 

: dip 



(13) 



(14) 



in which P k = P k (x, s) and r is the minimal time needed for reaching a fixed point P£ = PL*(x, s, r). Integral ( fl4| ) 
with velocity (13) yields 



P^x,s,r) 



v akT 



-l/u 



(15) 



where v = v k {s) 



°s+ao ■ Taking the algebraic transform inverse to (1 1 ), we find 



-1/1/ 



(16) 



Exponential renormalization [24, 26] corresponds to the limit s — > oo, at which linis^oo Vk{s) = 0, limg^oo sv k {s) 
a k — ao- Then (ffj ) gives 



lim pI(x,s,t) =p k ^i(x)exp ( — rx Qfe a ° 

ao 



(17) 



Accomplishing exponential renormalization of all sums appearing in expression of type (|17[), we follow the bootstrap 
procedure [|4|] according to the scheme p k {x) — > p* k {x, s, r) — > F k (x, t\,T2, r k ), with k > 1. 

Let us mention a recent innovation [^] that may improve significantly the convergence of the method, based on 
the determination of the control parameters from the knowledge of some moments of the function to reconstruct. Let 
us assume that we can obtain the first j — 1 moments /ij, i = 1, 2.., j of the sought function (f)(t),m some interval T, 



Mi 



f- L (f)(t)dt, 



(18) 
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so that for j=2 both zero and first moments are available etc... One can condition the control parameters T1T2 Tj as 
follows 

f*(t,T 1 T 2t ...,T J _ 1 )f- 1 dt = (19) 





Based on these conditions, one can attempt to solve two different problems, the first one corresponds to an approximate 
reconstruction of the function 4>(t) within the same interval [0,T] where moments are given or measured. The second 
problem consists in extrapolating to t > T. It is also possible to use an hybrid approach, where some controls are 
obtained from the agreement with expansion, while the remaining ones are found from the conditions on moments. 

4 Resummation of lower order expansions in o\. 

4.1 What can be extracted from the expansion of K^j 

Assume that the following extremely short expansion has been obtained, 

K{a Y )~l-aa Y , a =\~\ (4^0). (20) 

Consider the case a > 0. In order to find the behavior of K(ay) for arbitrary a Y , we continue it from the region of 

a Y — > self-similarly, along the most stable trajectory, with the crossover index s, determined by the condition of the 

minimum of the multiplier [EL E4l |29l] 

— — 1 + s 

m(ay , s) = 1 — aoy > (21) 



from where 

s(ay) = aoy (l — aoy) , ay < aT x l 2 , 



s — > 00, oy > a 1 ^ 2 , 



corresponding to the self-similar approximation 

s(ay) 



s(<7y) aa Y 

K *M = ( , "V' 2 ) = ( 2 " aa 2 Y )^~\ ay < a' 1 / 2 , (22) 



K*[a Y ) = exp(-acr y ), a Y > a' 1 ! 2 . (23) 
This suggests one self-similar expression (22) up to ayo = a -1 / 2 , and another (23), exponentially "soft", above this 



value. Strictly speaking, formula (22) is applicable for ay only up to (2/a) 1 / 2 , where it predicts a spurious zero of 
conductivity. In this particular case, the self-similar approximation plausibly reconstructs the exponential function for 
arbitrary ay, even in the absence of any a priori assumption on the asymptotic behavior at ay — > 00. 
For negative a (d > 2), there is no limiting value ayo and this yields 

siay ) \ s(aY) 



K*(ay)=i Tf' 2 , (24) 

\s{ay) + aa Y J 

which can be used for arbitrary ay. This expression appears to be more stable than the previously proposed exponen- 
tial solution exp(— aciy) for arbitrary ay, as indicated by the analysis of the multipliers. Expression (24) also predict 



a smaller conductivity than the exponential function for arbitrary ay suggesting that, for the most interesting case 
d = 3, the ansatz equation ( ||) should be replaced. Analysis of higher-order expansions will provide more details on 
the sought function. 
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4.2 Resummation of the second-order expansion in o Y . One-parameter formula 

Available from[|3|] in the next order in a Y , we have 



K 



[2] 
ef 



n ef + ^ ef 



K a 



1 + 



2 d Y 2\2 d Y 



(25) 



Below, for simplicity, we apply our resummation technique to the dimensionless quantity 



K(z) ~ 1 + a\z + ci2Z 2 , z = a Y , a\ 







,a 2 = i 








-a 




V2 





(26) 



Application of accuracy-through-order conditions, or of the superexponential approximants give, almost trivially, 
an exponential solution. Note that the Pade approximant available in this case, 



P(z) 



1 + (-02/01 + a\)z 



(27) 



1 — 02/01 z 

for d = 3, possesses a singularity at z = 12, which is wrong. 

The set of approximations to K(z), including the two starting terms from ( 26), can be written down as follows 

#o = l, 
Kx = 1 + aiz, 

and the expression for the renormalized quantity a\ can be readily obtained: 



s\ — a\z 



z Sl (z — > 00) 



(28) 



where the stabilizer 9\ should be negative, if we want to reproduce in the limit of z — ► 00, the correct, supposedly 
power-low behavior of the conductivity. A different set of approximations, which does not include the constant term 
from (26) into the renormalization procedure, has the form: 

K\ = a±z, 



K 2 = a\z + a 2 z 



(29) 
(30) 

(31) 



and applying the standard procedure of [|24|, g5J, we obtain 

#2=1 + a lZ [l - -7^]- (1+S2) => (-J-T- ■ 
ai(l + s 2 ) l + s 2 

Demanding now that both (28) and (|J) have the same power-law behavior at z — > 00, we find that 

S2 = S\ = S 

Requiring now the fulfillment of the stability criteria for the two available approximations in the form of the minimal- 
difference condition (see section ||), we obtain the condition that the negative stabilizer s should be determined from 
the minimum of the expression: 

r , \ -n+si , 

s 



\ -(1+*) 

l + sj 1 



(32) 



Generally speaking, it is sufficient to ask for an extremum of this difference. 



8 



In the case of d = 3, the maximum is located at the point s = —1.218. The final formulae have the following 
form: s 

Kl(z) = (—?—) , (33) 
V s — a\z I 



K 2 {z) = 1 + aiz 



1 a 2 z 



(34) 



ai(l + s) 

This last formula gives a lower bound for conductivity while the upper bound is simply exp{a\z) (equation (||)) which, 
in this case, is the only available "factor"-approximant based on all available (three) terms from the expansion. The 
corresponding multiplier is 

MHz) = + + ( ! _ -* x ) ~ iS+1) (35) 
2V ' oi(l + s)-022 V ai(l + a) J 

and the weighted average [ pT} , |9|] is given by 

c( Z) = i±jmmml_ i , (36) 

exp(-a 1 z) + |M 2 *(z)r 1 

providing the one-parametric formula for 3d-conductivity. See Fig.l for a comparison of different formulas for the 
conductivity as a function of the variance z = a\ defined in equation (^). The solid line corresponds to the average 
C(z), while the dotted line presents The dashed line is the celebrated exponential Landau-Lifshitz-Matheron 

(LLM) conjecture. The dash-dotted line corresponds to the result of resummation based on first-order expansion, 

K*(a Y ). 

4.3 Resummation of the third-order expansion in o\ . Two-parameters formula 

The expansion to the next order is given by [BJ|, 

K{z) ~ l + a lZ + a 2 z 2 + a 3 z 3 , a 3 {Z) = i Q - -Z, (37) 

where Z = 0.0042/3 (in the case of a Gaussian covariance ), or Z = 0.0014/3 (in the case of an exponentially 
decaying covariance). 

Using all terms from the expansion, we can create the following "odd" factor- approximant [[TT|], 

k*( r\ u (y \- {s2[Z)+1) (r . 4-as{Z)a L 

K 3 (z, Z) = 1 + a lZ 1 r^r~~iT ) ' S2 X =" 2 1 — o 1 7 \ > ^ 

V ai(s 2 (Z) + l) / a^-2a 3 (Z)ai 

while K^{z) = exp(aiz), which recovers expression (|8|). The approximant K 2 (z) gives an upper bound for the 
conductivity coefficient, while K£(z, Z) given by (38) provides a lower bound. The corresponding multipliers can be 
readily written down, 

_ A _ «V<"-<™ , (39) 

oi(l + s 2 (Z)) - a 2 z V a 1 (s 2 (Z) + l) ) 

M 2 {z) =exp(aiz) . (40) 
This allows us to obtain the weighted average of the conductivity coefficient 

= l + if^^lM^Z)!- 1 
V ' 7 exp(-oiz) + |M 3 *(z,Z)r 1 ' 
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providing a two-parameters formula for the 3d-conductivity. The results for C(z, Z) for Gaussian and exponential 
covariances are shown in Fig. 2, with the dotted line for the exponential case and with the solid line for the Gaussian 
case. These results are compared with the Landau-Lifshitz-Matheron (LLM) shown with dashed line. 



Numerical data in the exponential case are available till z = 7 Qlq ] and in the Gaussian case up to z = 6 [|5|]. Our 
results suggest that all formulas based on the expansion on a\ — up to the third order underestimate the conductivity 
and more terms are needed to improve the accuracy of the resummed expressions. 

A different approach aimed at increasing the accuracy consists in getting expressions for the conductivity in the 
limit of Oy — > oo, for instance from expansions in inverse powers of Oy, It is known (see Ref.[j8, 28, ^9|]) that when 
the asymptotic form of the solution is known, even only qualitatively, the formulas for the sought function can be 
improved very significantly. Even the knowledge of the leading power in the limit of large Oy would be of utmost 
importance. 

5 1/ d-expansion and resummation 
5.1 One-parametric case (d=3) 

Expression (^6|) for K(z) in the second order of perturbation theory can be re-written in the form of an expansion in 
the parameter 1/d with coefficients dependent on z, 

K(m) ~b (z) + b 1 (z)m + b 2 (z)m 2 , m=l/d; (42) 
b,{z) = l+ Z -+ Z ^ bx(z) = -z-^, &2(*) = y- (43) 



The theory of self-similar super-exponential approximants [ 24 , E6L P5L NOp then provides the following approximant 



h fb 2 \\ b\ 

—nmexp —Tim , n = 1, ti = 1 — — — 



K 2 {m) = b exp I — nmexp y—T 2 mJ J , n = 1, r 2 = 1 - . (44) 

iv"|(m) is located within the bounds given by K%(z) and exp(aiz) and provides a one-parametric formula for the 
d = 3-conductivity, as shown in Fig. 3. K% (m) (dotted line) appears to be located within the bounds outlined by C{z) 
(dashed) and exp(a\z) (solid line) and provides a one-parametric formula for the conductivity in three dimensions. 
The perturbative expression K(m) (dashed-dot line) is shown as well for comparison. 

5.2 Two-parametric case (d=3) 



Expression for K(z) (37) in the third order of perturbation theory can also be re-written in the form of an 1/d- 



expansion with coefficients dependent on z and Z, 

K(m, Z) ~ b (z, Z) + b x {z)m + b 2 (z)m 2 + b 3 (z)m 3 ; (45) 

2 2 3 2 3 

b (z,Z) = l + ^+ Z - + (±-Z)z 3 , bl {z) = -z- Z -- Z -, b 2 (z) = Z Y + Z T , (46) 

sy 2 3 *y 3 

6a(*) = Y + j, Hz) = -j. (47) 

We apply the technique of self-similar superexponential function in its variant detailed in [|9|, [l0|] , giving the following 
approximants 

Ki,(m,Z,Ti,T 2 ) = 6 exp (^-nmexp ^r 2 m^ , (48) 

Ks(m,Z,Ti,T 2 ,T 3 ) = 6 exp f^mexp ^r 2 mexp (^ rsm ))) ' ( - 49 - ) 
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1 (-bj + 6b 2 hb 3 - 6r 2 b b f b 2 - 3Tlb 2 b 2 2 ) 
6 T 2 b 2 ] b 1 b- i 



T 3 = \2T~i ' ( 5 °) 



In order to check whether the sequence of K* converges, we study their mapping multipliers, M*(t, ri, r 2 , . . . Tj) 
defined as 

M*(m, Z, n, . . . , Tj) = 1 ^ K*(m, Z,r u ..., Tj ) . (51) 

This definition of the multipliers allows us to compare the convergence of the expansion and of the renormalized 
expressions, making clear what can be expected a priori. 

This provides a matrix of self-similar approximants, indexed by the order j and by the number of control param- 
eters, 

K* 2l (m, Z) = K* (m, Z, n ,l), K* 22 (m, Z) = K* (m, Z,t x ,t 2 ), (52) 

K* 1 (m,Z)=K* 3 (m,Z, n , 1,1), K* 32 (m, Z) = K* 3 (m, Z, n , r 2 , 1), K* 33 (m, Z) = K* 3 (m, Z, n ,r 2 ,r 3 ), (53) 
Approximants K 22 (m, Z) and K 32 (m, Z) form the closest pair. Their average 

K * {m z) = gjgK Z) \M* 2 (m, ZT 1 + K* 2 (m, Z) \M* 2 (m, Z)^ 
' " \M* 2 (m,Z)\- l + \M* 2 (m,Z)\- 1 

is located within the bounds given by K 3 {z, Z) ([38|) and exp(aiz) and provides a useful formula for the conductivity 
coefficient. The results for the exponential and Gaussian covariances are shown in Figures 4 and 5 respectively. 
In the exponential case, good agreement between K*(m,Z) (dash-dot) and the Landau -Lifshitz-Matheron (LLM) 
conjecture (solid) remains valid till z w 11. The average behavior appears to be located within the bounds outlined 
by K 22 (m, Z) (dashed) and K 32 (m, Z) (dotted) and provides a reasonable formula for the conductivity. Numerical 



data are available till z = 7 Q16| ] and they agree well with our lower bound. 

In the Gaussian case K 32 (m, Z) (dotted line) gives an approximation which is closer to the Landau-Lifshitz- 
Matheron (LLM) conjecture (solid) than K 22 (m, Z) (dashed line). In this case, numerical results are available up to 
* = 6|§. 

We conclude, tentatively, that all formulas based on 1/d— expansion in the third order provide rather accurate 
expressions for conductivity for small and moderate variances and disagree with the LLM-conjecture for very large 
variances. To the best of our knowledge, the region of large variances is not accessible by other techniques, numer- 
ically or theoretically. All formulas based on the novel proposed 1 / d— expansion in the third order provide rather 
accurate expressions for the conductivity coefficient. We note that 1 / d-expansions may be faster converging than the 
original expansion in variances and should be investigated future. 



6 Future directions 

We have shown that it is possible to extend the moment equation approach using the self-similar functional renormal- 
ization method in order to provide a stable and robust estimation of the transport variables of the three-dimensional 
medium at large values of the perturbation parameter <Ty. 

Future directions include the extension of the self-similar functional renormalization method to go beyond the 
characterization of the heterogeneity solely in terms of the variance and consider also the dependence of the transport 
properties with respect to the skewness (third normalized cumulant) and kurtosis (fourth normalized cumulant). 

Acknowledgment: We are grateful to D.M. Tartakovsky for introducing us into the subject and for useful discus- 
sions. 
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Figures Captions 



Figure 1. Dependence of different estimators of the conductivity as a function of the variance z = a Y defined 
in equation (|26|). The solid line corresponds to C(z), the dotted line shows K^z). The dashed line is the Landau- 
Lifshitz-Matheron (LLM) conjecture. K*(cjy) (dash-dotted line) corresponds to the result of resummation based on 
the first-order expansion. 



Figure 2. Dependence of the conductivity C(z, Z) as a function of the variance z = a Y defined in equation 
(26) for Gaussian and exponential covariances, shown with solid and dotted lines respectively. For comparison, the 
Landau-Lifshitz-Matheron (LLM) conjecture is shown in dashed lines. 



Figure 3. As a function of the variance z = a Y defined in equation (|^), this figure shows a comparison between 
the conductivity if|(m) (dotted line) with C{z) (dashed) and with exp(aiz) (solid line). The Perturbative expansion 
K{m) is shown with the dashed-dot line. 



Figure 4. Exponential covariance: K*(m, Z) (dash-dot) is compared with the Landau-Lifshitz-Matheron (LLM) 
conjecture (solid line). Two approximants K% 2 (m, Z) (dashed) and K^ 2 {m, Z) (dotted) are shown as well. 



Figure 5. Gaussian case: Approximant iv"| 2 (m, Z) (dotted) compared with the Landau-Lifshitz-Matheron (LLM) 
conjecture (solid line). The approximant K^im, Z) (dashed line) is shown as well. 
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